#library(R2jags)
library(R2jags)
library(tidyverse)
library(lubridate)
library(knitr)
library(broom)
leaders_data <- get(load("Leaders.RData"))
leaders_data$Censored <- as.factor(leaders_data$Censored)
leaders_data$Type <- as.factor(leaders_data$Type)
leaders_data$Age.Event <- as.numeric(leaders_data$Age.Event)
summary(leaders_data)
## Name Birth.Date Event.Date Age.Event
## Length:177 Min. :1324-05-14 Min. :1368-03-29 Min. : 9.259
## Class :character 1st Qu.:1507-09-16 1st Qu.:1572-07-05 1st Qu.:53.985
## Mode :character Median :1705-10-31 Median :1758-05-03 Median :67.529
## Mean :1674-01-23 Mean :1737-03-27 Mean :63.170
## 3rd Qu.:1833-08-20 3rd Qu.:1886-11-18 3rd Qu.:78.242
## Max. :1961-08-04 Max. :2020-07-31 Max. :95.830
## Censored Type Fail centdate time
## 0:166 ChinaEmp :26 Min. :0.0000 Min. :1300-01-01 Min. :0.2436
## 1: 11 DalaiLama:14 1st Qu.:1.0000 1st Qu.:1300-01-01 1st Qu.:2.0770
## JapanEmp :30 Median :1.0000 Median :1300-01-01 Median :4.0582
## Pope :63 Mean :0.9379 Mean :1300-01-01 Mean :3.7406
## UsPres :44 3rd Qu.:1.0000 3rd Qu.:1300-01-01 3rd Qu.:5.3362
## Max. :1.0000 Max. :1300-01-01 Max. :6.6157
ggplot(data = leaders_data) +
geom_boxplot(aes(x = Type, y = Age.Event))
ggplot(data = leaders_data) +
geom_point(aes(x = Birth.Date, y = Age.Event))
\[ T_i \sim Weibull(\alpha, \mu) \\ log(\mu) = \beta_0 + \beta_1x_{i1}+\beta_2x_{i2}+...+\beta_px_{ip}\\ \] Log of mu is done so that the mu parameter is always positive (the support on mu in Weibull is that mu > 0)
\[ T_i \sim Weibull(r, \mu_i) \\ \log(\mu_i) = \beta_0+\beta_1 \text{(Year of Birth_i)}+\beta_2 \text{(Leadership_i = UsPres)} + \beta_3\text{(Leadership_i = ChinaEmp)} \\ + \beta_4\text{(Leadership_i = DalaiLama)} +\beta_5\text{(Leadership_i = JapanEmp)} \\ + \beta_6\text{(Year of Birth_i * (Leadership_i = UsPres))} \\ + \beta_7\text{(Year of Birth_i * (Leadership_i = ChinaEmp))} \\ + \beta_8\text{(Year of Birth_i * (Leadership_i = DalaiLama))} \\ + \beta_9\text{(Year of Birth_i * (Leadership_i = JapanEmp))} \]
leaders_data$Birth.Year <- year(leaders_data$Birth.Date)
leaders_data <- leaders_data %>%
mutate(TypeChinaEmp = as.factor(case_when(Type == "ChinaEmp" ~ 1,
TRUE ~ 0)),
TypeDalaiLama = as.factor(case_when(Type == "DalaiLama" ~ 1,
TRUE ~ 0)),
TypeJapanEmp = as.factor(case_when(Type == "JapanEmp" ~ 1,
TRUE ~ 0)),
TypePope = as.factor(case_when(Type == "Pope" ~ 1,
TRUE ~ 0)),
TypeUsPres = as.factor(case_when(Type == "UsPres" ~ 1,
TRUE ~ 0)))
leaders_nopred <- leaders_data %>%
filter(!(Name %in% c("14th Dalai Lama (Tenzin Gyatso)", "Naruhito (Emperor Reiwa)", "Barack Obama", "Francis")))
survival_raw <- leaders_nopred$Age.Event
n <- length(survival_raw)
censored <- leaders_nopred$Censored
survival <- ifelse(censored == 0, survival_raw, NA)
censoring_limits <- ifelse(censored == 0, 120, survival_raw)
x_1 <- leaders_nopred$Birth.Year - mean(leaders_nopred$Birth.Year)
x_2 <- leaders_nopred$TypeUsPres
x_3 <- leaders_nopred$TypeChinaEmp
x_4 <- leaders_nopred$TypeDalaiLama
x_5 <- leaders_nopred$TypeJapanEmp
# add interaction b/w Birth.Year and Type
x_6 <- x_1*as.numeric(as.character(x_2))
x_7 <- x_1*as.numeric(as.character(x_3))
x_8 <- x_1*as.numeric(as.character(x_4))
x_9 <- x_1*as.numeric(as.character(x_5))
set.seed(123)
second_approach <- function() {
for(i in 1:n_censored) {
z_censored[i] ~ dpois(phi_censored[i])
phi_censored[i] <- mu_censored[i] * pow(t_censored[i], r)
mu_censored[i] <- exp(beta_censored[i])
beta_censored[i] <- beta_0 + beta_1*x_1_censored[i] + beta_2*x_2_censored[i] + beta_3*x_3_censored[i] + beta_4*x_4_censored[i] + beta_5*x_5_censored[i] + beta_6*x_6_censored[i] + beta_7*x_7_censored[i] + beta_8*x_8_censored[i] + beta_9*x_9_censored[i]
}
for(j in 1:n_non_censored) { #total rows - 7 ** 7 are censored, rest are not censored
survival_non_censored[j] ~ dweib(r, mu[j])
mu[j] <- exp(beta[j])
beta[j] <- beta_0 + beta_1*x_1_non_censored[j] + beta_2*x_2_non_censored[j] + beta_3*x_3_non_censored[j] + beta_4*x_4_non_censored[j] + beta_5*x_5_non_censored[j] + beta_6*x_6_non_censored[j] + beta_7*x_7_non_censored[j] + beta_8*x_8_non_censored[j] + beta_9*x_9_non_censored[j]
}
beta_0 ~ dnorm(0.0, 1.0E-3) # Prior on beta_0 is normal with low precision
beta_1 ~ dnorm(0.0, 1.0E-3) # Prior on beta_1 is normal with low precision
beta_2 ~ dnorm(0.0, 1.0E-3) # Prior on beta_2 is normal with low precision
beta_3 ~ dnorm(0.0, 1.0E-3) # Prior on beta_3 is normal with low precision
beta_4 ~ dnorm(0.0, 1.0E-3) # Prior on beta_4 is normal with low precision
beta_5 ~ dnorm(0.0, 1.0E-3) # Prior on beta_5 is normal with low precision
beta_6 ~ dnorm(0.0, 1.0E-3) # Prior on beta_6 is normal with low precision
beta_7 ~ dnorm(0.0, 1.0E-3) # Prior on beta_7 is normal with low precision
beta_8 ~ dnorm(0.0, 1.0E-3) # Prior on beta_8 is normal with low precision
beta_9 ~ dnorm(0.0, 1.0E-3) # Prior on beta_9 is normal with low precision
r ~ dexp(1) # Prior on r
# Predictive distribution of age at the new values
beta_Francis <- beta_0 + beta_1*year_Francis
mu_Francis <- exp(beta_Francis)
survival_Francis ~ dweib(r, mu_Francis) %_% T(present_length_Francis, upper_length)
age_Francis_predictive <- survival_Francis
beta_Obama <- beta_0 + (beta_1+beta_6)*year_Obama + beta_2
mu_Obama <- exp(beta_Obama)
survival_Obama ~ dweib(r, mu_Obama) %_% T(present_length_Obama, upper_length)
age_Obama_predictive <- survival_Obama
beta_Dalai <- beta_0 + (beta_1+beta_8)*year_Dalai + beta_4
mu_Dalai <- exp(beta_Dalai)
survival_Dalai ~ dweib(r, mu_Dalai) %_% T(present_length_Dalai, upper_length)
age_Dalai_predictive <- survival_Dalai
beta_Naruhito <- beta_0 + (beta_1+beta_9)*year_Naruhito + beta_5
mu_Naruhito <- exp(beta_Naruhito)
survival_Naruhito ~ dweib(r, mu_Naruhito) %_% T(present_length_Naruhito, upper_length)
age_Naruhito_predictive <- survival_Naruhito
}
censored_dex = which(censored==1) # index of ppl who are censored
z_censored <- rep(0, length(censored_dex))
t_censored <- censoring_limits[censored_dex]
x_1_censored <- x_1[censored_dex]
x_2_censored <- x_2[censored_dex]
x_3_censored <- x_3[censored_dex]
x_4_censored <- x_4[censored_dex]
x_5_censored <- x_5[censored_dex]
x_6_censored <- x_6[censored_dex]
x_7_censored <- x_7[censored_dex]
x_8_censored <- x_8[censored_dex]
x_9_censored <- x_9[censored_dex]
n_censored <- length(censored_dex)
survival_non_censored <- survival[-censored_dex] # Remove ALL ALIVE PEOPLE
x_1_non_censored <- x_1[-censored_dex]
x_2_non_censored <- x_2[-censored_dex]
x_3_non_censored <- x_3[-censored_dex]
x_4_non_censored <- x_4[-censored_dex]
x_5_non_censored <- x_5[-censored_dex]
x_6_non_censored <- x_6[-censored_dex]
x_7_non_censored <- x_7[-censored_dex]
x_8_non_censored <- x_8[-censored_dex]
x_9_non_censored <- x_9[-censored_dex]
n_non_censored <- length(survival_non_censored)
birth_year_Francis <- leaders_data$Birth.Year[leaders_data$Name == "Francis"]
year_Francis <- birth_year_Francis - mean(leaders_nopred$Birth.Year)
present_length_Francis <- leaders_data$Age.Event[leaders_data$Name =="Francis"]
birth_year_Obama <- leaders_data$Birth.Year[leaders_data$Name == "Barack Obama"]
year_Obama <- birth_year_Obama - mean(leaders_nopred$Birth.Year)
present_length_Obama <- leaders_data$Age.Event[leaders_data$Name =="Barack Obama"]
birth_year_Dalai <- leaders_data$Birth.Year[leaders_data$Name == "14th Dalai Lama (Tenzin Gyatso)"]
year_Dalai <- birth_year_Dalai - mean(leaders_nopred$Birth.Year)
present_length_Dalai <- leaders_data$Age.Event[leaders_data$Name =="14th Dalai Lama (Tenzin Gyatso)"]
birth_year_Naruhito <- leaders_data$Birth.Year[leaders_data$Name == "Naruhito (Emperor Reiwa)"]
year_Naruhito <- birth_year_Naruhito - mean(leaders_nopred$Birth.Year)
present_length_Naruhito <- leaders_data$Age.Event[leaders_data$Name =="Naruhito (Emperor Reiwa)"]
upper_length <- 120
data_Popes_second_approach <- list("n_censored",
"n_non_censored",
"z_censored",
"t_censored",
"x_1_censored",
"x_2_censored",
"x_3_censored",
"x_4_censored",
"x_5_censored",
"x_6_censored",
"x_7_censored",
"x_8_censored",
"x_9_censored",
"survival_non_censored",
"x_1_non_censored",
"x_2_non_censored",
"x_3_non_censored",
"x_4_non_censored",
"x_5_non_censored",
"x_6_non_censored",
"x_7_non_censored",
"x_8_non_censored",
"x_9_non_censored",
"year_Francis",
"year_Obama",
"year_Dalai",
"year_Naruhito",
"present_length_Francis",
"present_length_Obama",
"present_length_Dalai",
"present_length_Naruhito",
"upper_length")
Bayesian_Popes_second_approach <- jags(data = data_Popes_second_approach,
parameters.to.save = c("beta_0",
"beta_1",
"beta_2",
"beta_3",
"beta_4",
"beta_5",
"beta_6",
"beta_7",
"beta_8",
"beta_9",
"r",
"age_Francis_predictive",
"age_Obama_predictive",
"age_Dalai_predictive",
"age_Naruhito_predictive",
"survival_non_censored"
),
n.iter = 50000,
n.chains = 3,
model.file = second_approach)
## module glm loaded
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 173
## Unobserved stochastic nodes: 15
## Total graph size: 2385
##
## Initializing model
Bayesian_Popes_second_approach
## Inference for Bugs model at "/tmp/RtmpyN0MZP/model49d3f84c442.txt", fit using jags,
## 3 chains, each with 50000 iterations (first 25000 discarded), n.thin = 25
## n.sims = 3000 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75%
## age_Dalai_predictive 100.942 9.940 85.743 92.172 100.329 109.281
## age_Francis_predictive 104.921 10.069 85.599 96.968 106.165 113.568
## age_Naruhito_predictive 98.891 15.455 66.047 87.555 101.833 112.140
## age_Obama_predictive 99.278 15.361 64.529 89.100 102.630 112.128
## beta_0 -22.776 1.596 -25.941 -23.801 -22.745 -21.742
## beta_1 -0.002 0.001 -0.003 -0.002 -0.002 -0.001
## beta_2 0.601 0.440 -0.282 0.307 0.610 0.906
## beta_3 1.324 0.268 0.790 1.148 1.324 1.509
## beta_4 1.694 0.335 1.019 1.474 1.707 1.926
## beta_5 0.716 0.244 0.241 0.556 0.712 0.881
## beta_6 -0.002 0.002 -0.007 -0.004 -0.002 0.000
## beta_7 0.001 0.001 -0.002 0.000 0.001 0.002
## beta_8 0.007 0.002 0.004 0.006 0.007 0.008
## beta_9 0.000 0.001 -0.002 -0.001 0.000 0.001
## r 4.240 0.265 3.743 4.045 4.241 4.419
## survival_non_censored[1] 84.868 0.000 84.868 84.868 84.868 84.868
## survival_non_censored[2] 65.947 0.000 65.947 65.947 65.947 65.947
## survival_non_censored[3] 80.857 0.000 80.857 80.857 80.857 80.857
## survival_non_censored[4] 81.517 0.000 81.517 81.517 81.517 81.517
## survival_non_censored[5] 82.601 0.000 82.601 82.601 82.601 82.601
## survival_non_censored[6] 81.695 0.000 81.695 81.695 81.695 81.695
## survival_non_censored[7] 67.168 0.000 67.168 67.168 67.168 67.168
## survival_non_censored[8] 79.214 0.000 79.214 79.214 79.214 79.214
## survival_non_censored[9] 93.380 0.000 93.380 93.380 93.380 93.380
## survival_non_censored[10] 85.736 0.000 85.736 85.736 85.736 85.736
## survival_non_censored[11] 80.698 0.000 80.698 80.698 80.698 80.698
## survival_non_censored[12] 69.027 0.000 69.027 69.027 69.027 69.027
## survival_non_censored[13] 68.468 0.000 68.468 68.468 68.468 68.468
## survival_non_censored[14] 81.013 0.000 81.013 81.013 81.013 81.013
## survival_non_censored[15] 81.676 0.000 81.676 81.676 81.676 81.676
## survival_non_censored[16] 68.893 0.000 68.893 68.893 68.893 68.893
## survival_non_censored[17] 75.907 0.000 75.907 75.907 75.907 75.907
## survival_non_censored[18] 83.088 0.000 83.088 83.088 83.088 83.088
## survival_non_censored[19] 87.830 0.000 87.830 87.830 87.830 87.830
## survival_non_censored[20] 81.049 0.000 81.049 81.049 81.049 81.049
## survival_non_censored[21] 68.816 0.000 68.816 68.816 68.816 68.816
## survival_non_censored[22] 71.652 0.000 71.652 71.652 71.652 71.652
## survival_non_censored[23] 85.541 0.000 85.541 85.541 85.541 85.541
## survival_non_censored[24] 80.780 0.000 80.780 80.780 80.780 80.780
## survival_non_censored[25] 78.242 0.000 78.242 78.242 78.242 78.242
## survival_non_censored[26] 86.026 0.000 86.026 86.026 86.026 86.026
## survival_non_censored[27] 69.864 0.000 69.864 69.864 69.864 69.864
## survival_non_censored[28] 68.268 0.000 68.268 68.268 68.268 68.268
## survival_non_censored[29] 80.674 0.000 80.674 80.674 80.674 80.674
## survival_non_censored[30] 76.315 0.000 76.315 76.315 76.315 76.315
## survival_non_censored[31] 69.478 0.000 69.478 69.478 69.478 69.478
## survival_non_censored[32] 70.366 0.000 70.366 70.366 70.366 70.366
## survival_non_censored[33] 69.862 0.000 69.862 69.862 69.862 69.862
## survival_non_censored[34] 69.021 0.000 69.021 69.021 69.021 69.021
## survival_non_censored[35] 72.446 0.000 72.446 72.446 72.446 72.446
## survival_non_censored[36] 56.676 0.000 56.676 56.676 56.676 56.676
## survival_non_censored[37] 69.147 0.000 69.147 69.147 69.147 69.147
## survival_non_censored[38] 68.704 0.000 68.704 68.704 68.704 68.704
## survival_non_censored[39] 83.255 0.000 83.255 83.255 83.255 83.255
## survival_non_censored[40] 68.287 0.000 68.287 68.287 68.287 68.287
## survival_non_censored[41] 66.691 0.000 66.691 66.691 66.691 66.691
## survival_non_censored[42] 83.135 0.000 83.135 83.135 83.135 83.135
## survival_non_censored[43] 53.985 0.000 53.985 53.985 53.985 53.985
## survival_non_censored[44] 67.529 0.000 67.529 67.529 67.529 67.529
## survival_non_censored[45] 81.695 0.000 81.695 81.695 81.695 81.695
## survival_non_censored[46] 56.331 0.000 56.331 56.331 56.331 56.331
## survival_non_censored[47] 64.534 0.000 64.534 64.534 64.534 64.534
## survival_non_censored[48] 45.971 0.000 45.971 45.971 45.971 45.971
## survival_non_censored[49] 69.213 0.000 69.213 69.213 69.213 69.213
## survival_non_censored[50] 64.386 0.000 64.386 64.386 64.386 64.386
## survival_non_censored[51] 72.624 0.000 72.624 72.624 72.624 72.624
## survival_non_censored[52] 60.564 0.000 60.564 60.564 60.564 60.564
## survival_non_censored[53] 70.062 0.000 70.062 70.062 70.062 70.062
## survival_non_censored[54] 54.418 0.000 54.418 54.418 54.418 54.418
## survival_non_censored[55] 58.823 0.000 58.823 58.823 58.823 58.823
## survival_non_censored[56] 79.595 0.000 79.595 79.595 79.595 79.595
## survival_non_censored[57] 57.355 0.000 57.355 57.355 57.355 57.355
## survival_non_censored[58] 64.142 0.000 64.142 64.142 64.142 64.142
## survival_non_censored[59] 62.133 0.000 62.133 62.133 62.133 62.133
## survival_non_censored[60] 91.135 0.000 91.135 91.135 91.135 91.135
## survival_non_censored[61] 67.844 0.000 67.844 67.844 67.844 67.844
## survival_non_censored[62] 67.808 0.000 67.808 67.808 67.808 67.808
## survival_non_censored[63] 90.675 0.000 90.675 90.675 90.675 90.675
## survival_non_censored[64] 83.222 0.000 83.222 83.222 83.222 83.222
## survival_non_censored[65] 85.284 0.000 85.284 85.284 85.284 85.284
## survival_non_censored[66] 73.180 0.000 73.180 73.180 73.180 73.180
## survival_non_censored[67] 80.619 0.000 80.619 80.619 80.619 80.619
## survival_non_censored[68] 78.231 0.000 78.231 78.231 78.231 78.231
## survival_non_censored[69] 79.630 0.000 79.630 79.630 79.630 79.630
## survival_non_censored[70] 68.145 0.000 68.145 68.145 68.145 68.145
## survival_non_censored[71] 71.806 0.000 71.806 71.806 71.806 71.806
## survival_non_censored[72] 53.615 0.000 53.615 53.615 53.615 53.615
## survival_non_censored[73] 65.618 0.000 65.618 65.618 65.618 65.618
## survival_non_censored[74] 74.163 0.000 74.163 74.163 74.163 74.163
## survival_non_censored[75] 64.873 0.000 64.873 64.873 64.873 64.873
## survival_non_censored[76] 77.106 0.000 77.106 77.106 77.106 77.106
## survival_non_censored[77] 56.170 0.000 56.170 56.170 56.170 56.170
## survival_non_censored[78] 66.585 0.000 66.585 66.585 66.585 66.585
## survival_non_censored[79] 63.239 0.000 63.239 63.239 63.239 63.239
## survival_non_censored[80] 70.289 0.000 70.289 70.289 70.289 70.289
## survival_non_censored[81] 49.834 0.000 49.834 49.834 49.834 49.834
## survival_non_censored[82] 57.120 0.000 57.120 57.120 57.120 57.120
## survival_non_censored[83] 71.266 0.000 71.266 71.266 71.266 71.266
## survival_non_censored[84] 67.559 0.000 67.559 67.559 67.559 67.559
## survival_non_censored[85] 58.623 0.000 58.623 58.623 58.623 58.623
## survival_non_censored[86] 60.192 0.000 60.192 60.192 60.192 60.192
## survival_non_censored[87] 72.474 0.000 72.474 72.474 72.474 72.474
## survival_non_censored[88] 67.097 0.000 67.097 67.097 67.097 67.097
## survival_non_censored[89] 57.744 0.000 57.744 57.744 57.744 57.744
## survival_non_censored[90] 60.504 0.000 60.504 60.504 60.504 60.504
## survival_non_censored[91] 90.193 0.000 90.193 90.193 90.193 90.193
## survival_non_censored[92] 63.195 0.000 63.195 63.195 63.195 63.195
## survival_non_censored[93] 88.632 0.000 88.632 88.632 88.632 88.632
## survival_non_censored[94] 78.450 0.000 78.450 78.450 78.450 78.450
## survival_non_censored[95] 46.483 0.000 46.483 46.483 46.483 46.483
## survival_non_censored[96] 64.405 0.000 64.405 64.405 64.405 64.405
## survival_non_censored[97] 81.227 0.000 81.227 81.227 81.227 81.227
## survival_non_censored[98] 93.451 0.000 93.451 93.451 93.451 93.451
## survival_non_censored[99] 93.328 0.000 93.328 93.328 93.328 93.328
## survival_non_censored[100] 94.467 0.000 94.467 94.467 94.467 94.467
## survival_non_censored[101] 69.673 0.000 69.673 69.673 69.673 69.673
## survival_non_censored[102] 64.277 0.000 64.277 64.277 64.277 64.277
## survival_non_censored[103] 46.782 0.000 46.782 46.782 46.782 46.782
## survival_non_censored[104] 35.877 0.000 35.877 35.877 35.877 35.877
## survival_non_censored[105] 36.235 0.000 36.235 36.235 36.235 36.235
## survival_non_censored[106] 28.476 0.000 28.476 28.476 28.476 28.476
## survival_non_censored[107] 39.751 0.000 39.751 39.751 39.751 39.751
## survival_non_censored[108] 34.856 0.000 34.856 34.856 34.856 34.856
## survival_non_censored[109] 29.478 0.000 29.478 29.478 29.478 29.478
## survival_non_censored[110] 60.186 0.000 60.186 60.186 60.186 60.186
## survival_non_censored[111] 35.337 0.000 35.337 35.337 35.337 35.337
## survival_non_censored[112] 56.956 0.000 56.956 56.956 56.956 56.956
## survival_non_censored[113] 38.081 0.000 38.081 38.081 38.081 38.081
## survival_non_censored[114] 21.769 0.000 21.769 21.769 21.769 21.769
## survival_non_censored[115] 33.243 0.000 33.243 33.243 33.243 33.243
## survival_non_censored[116] 22.897 0.000 22.897 22.897 22.897 22.897
## survival_non_censored[117] 68.624 0.000 68.624 68.624 68.624 68.624
## survival_non_censored[118] 56.816 0.000 56.816 56.816 56.816 56.816
## survival_non_censored[119] 87.203 0.000 87.203 87.203 87.203 87.203
## survival_non_censored[120] 59.800 0.000 59.800 59.800 59.800 59.800
## survival_non_censored[121] 67.444 0.000 67.444 67.444 67.444 67.444
## survival_non_censored[122] 30.100 0.000 30.100 30.100 30.100 30.100
## survival_non_censored[123] 18.710 0.000 18.710 18.710 18.710 18.710
## survival_non_censored[124] 37.251 0.000 37.251 37.251 37.251 37.251
## survival_non_censored[125] 61.689 0.000 61.689 61.689 61.689 61.689
## survival_non_censored[126] 82.998 0.000 82.998 82.998 82.998 82.998
## survival_non_censored[127] 66.998 0.000 66.998 66.998 66.998 66.998
## survival_non_censored[128] 44.999 0.000 44.999 44.999 44.999 44.999
## survival_non_censored[129] 27.907 0.000 27.907 27.907 27.907 27.907
## survival_non_censored[130] 64.999 0.000 64.999 64.999 64.999 64.999
## survival_non_censored[131] 23.707 0.000 23.707 23.707 23.707 23.707
## survival_non_censored[132] 49.002 0.000 49.002 49.002 49.002 49.002
## survival_non_censored[133] 45.996 0.000 45.996 45.996 45.996 45.996
## survival_non_censored[134] 9.259 0.000 9.259 9.259 9.259 9.259
## survival_non_censored[135] 21.506 0.000 21.506 21.506 21.506 21.506
## survival_non_censored[136] 17.248 0.000 17.248 17.248 17.248 17.248
## survival_non_censored[137] 18.242 0.000 18.242 18.242 18.242 18.242
## survival_non_censored[138] 57.843 0.000 57.843 57.843 57.843 57.843
## survival_non_censored[139] 40.241 0.000 40.241 40.241 40.241 40.241
## survival_non_censored[140] 51.652 0.000 51.652 51.652 51.652 51.652
## survival_non_censored[141] 77.352 0.000 77.352 77.352 77.352 77.352
## survival_non_censored[142] 56.331 0.000 56.331 56.331 56.331 56.331
## survival_non_censored[143] 27.302 0.000 27.302 27.302 27.302 27.302
## survival_non_censored[144] 51.526 0.000 51.526 51.526 51.526 51.526
## survival_non_censored[145] 58.300 0.000 58.300 58.300 58.300 58.300
## survival_non_censored[146] 63.493 0.000 63.493 63.493 63.493 63.493
## survival_non_censored[147] 62.667 0.000 62.667 62.667 62.667 62.667
## survival_non_censored[148] 75.639 0.000 75.639 75.639 75.639 75.639
## survival_non_censored[149] 45.818 0.000 45.818 45.818 45.818 45.818
## survival_non_censored[150] 84.203 0.000 84.203 84.203 84.203 84.203
## survival_non_censored[151] 72.903 0.000 72.903 72.903 72.903 72.903
## survival_non_censored[152] 21.528 0.000 21.528 21.528 21.528 21.528
## survival_non_censored[153] 47.220 0.000 47.220 47.220 47.220 47.220
## survival_non_censored[154] 78.209 0.000 78.209 78.209 78.209 78.209
## survival_non_censored[155] 34.237 0.000 34.237 34.237 34.237 34.237
## survival_non_censored[156] 35.318 0.000 35.318 35.318 35.318 35.318
## survival_non_censored[157] 30.300 0.000 30.300 30.300 30.300 30.300
## survival_non_censored[158] 21.380 0.000 21.380 21.380 21.380 21.380
## survival_non_censored[159] 73.248 0.000 73.248 73.248 73.248 73.248
## survival_non_censored[160] 21.363 0.000 21.363 21.363 21.363 21.363
## survival_non_censored[161] 69.216 0.000 69.216 69.216 69.216 69.216
## survival_non_censored[162] 45.936 0.000 45.936 45.936 45.936 45.936
## survival_non_censored[163] 35.526 0.000 35.526 35.526 35.526 35.526
## survival_non_censored[164] 59.734 0.000 59.734 59.734 59.734 59.734
## survival_non_censored[165] 47.316 0.000 47.316 47.316 47.316 47.316
## survival_non_censored[166] 87.693 0.000 87.693 87.693 87.693 87.693
## deviance 1423.309 4.725 1416.048 1419.895 1422.654 1426.055
## 97.5% Rhat n.eff
## age_Dalai_predictive 118.774 1.002 1100
## age_Francis_predictive 119.424 1.001 3000
## age_Naruhito_predictive 119.328 1.001 2900
## age_Obama_predictive 119.270 1.001 3000
## beta_0 -19.587 1.051 53
## beta_1 0.000 1.002 1100
## beta_2 1.416 1.013 160
## beta_3 1.848 1.009 250
## beta_4 2.331 1.022 110
## beta_5 1.201 1.012 170
## beta_6 0.003 1.003 740
## beta_7 0.003 1.001 2400
## beta_8 0.010 1.003 890
## beta_9 0.003 1.002 1700
## r 4.762 1.027 110
## survival_non_censored[1] 84.868 1.000 1
## survival_non_censored[2] 65.947 1.000 1
## survival_non_censored[3] 80.857 1.000 1
## survival_non_censored[4] 81.517 1.000 1
## survival_non_censored[5] 82.601 1.000 1
## survival_non_censored[6] 81.695 1.000 1
## survival_non_censored[7] 67.168 1.000 1
## survival_non_censored[8] 79.214 1.000 1
## survival_non_censored[9] 93.380 1.000 1
## survival_non_censored[10] 85.736 1.000 1
## survival_non_censored[11] 80.698 1.000 1
## survival_non_censored[12] 69.027 1.000 1
## survival_non_censored[13] 68.468 1.000 1
## survival_non_censored[14] 81.013 1.000 1
## survival_non_censored[15] 81.676 1.000 1
## survival_non_censored[16] 68.893 1.000 1
## survival_non_censored[17] 75.907 1.000 1
## survival_non_censored[18] 83.088 1.000 1
## survival_non_censored[19] 87.830 1.000 1
## survival_non_censored[20] 81.049 1.000 1
## survival_non_censored[21] 68.816 1.000 1
## survival_non_censored[22] 71.652 1.000 1
## survival_non_censored[23] 85.541 1.000 1
## survival_non_censored[24] 80.780 1.000 1
## survival_non_censored[25] 78.242 1.000 1
## survival_non_censored[26] 86.026 1.000 1
## survival_non_censored[27] 69.864 1.000 1
## survival_non_censored[28] 68.268 1.000 1
## survival_non_censored[29] 80.674 1.000 1
## survival_non_censored[30] 76.315 1.000 1
## survival_non_censored[31] 69.478 1.000 1
## survival_non_censored[32] 70.366 1.000 1
## survival_non_censored[33] 69.862 1.000 1
## survival_non_censored[34] 69.021 1.000 1
## survival_non_censored[35] 72.446 1.000 1
## survival_non_censored[36] 56.676 1.000 1
## survival_non_censored[37] 69.147 1.000 1
## survival_non_censored[38] 68.704 1.000 1
## survival_non_censored[39] 83.255 1.000 1
## survival_non_censored[40] 68.287 1.000 1
## survival_non_censored[41] 66.691 1.000 1
## survival_non_censored[42] 83.135 1.000 1
## survival_non_censored[43] 53.985 1.000 1
## survival_non_censored[44] 67.529 1.000 1
## survival_non_censored[45] 81.695 1.000 1
## survival_non_censored[46] 56.331 1.000 1
## survival_non_censored[47] 64.534 1.000 1
## survival_non_censored[48] 45.971 1.000 1
## survival_non_censored[49] 69.213 1.000 1
## survival_non_censored[50] 64.386 1.000 1
## survival_non_censored[51] 72.624 1.000 1
## survival_non_censored[52] 60.564 1.000 1
## survival_non_censored[53] 70.062 1.000 1
## survival_non_censored[54] 54.418 1.000 1
## survival_non_censored[55] 58.823 1.000 1
## survival_non_censored[56] 79.595 1.000 1
## survival_non_censored[57] 57.355 1.000 1
## survival_non_censored[58] 64.142 1.000 1
## survival_non_censored[59] 62.133 1.000 1
## survival_non_censored[60] 91.135 1.000 1
## survival_non_censored[61] 67.844 1.000 1
## survival_non_censored[62] 67.808 1.000 1
## survival_non_censored[63] 90.675 1.000 1
## survival_non_censored[64] 83.222 1.000 1
## survival_non_censored[65] 85.284 1.000 1
## survival_non_censored[66] 73.180 1.000 1
## survival_non_censored[67] 80.619 1.000 1
## survival_non_censored[68] 78.231 1.000 1
## survival_non_censored[69] 79.630 1.000 1
## survival_non_censored[70] 68.145 1.000 1
## survival_non_censored[71] 71.806 1.000 1
## survival_non_censored[72] 53.615 1.000 1
## survival_non_censored[73] 65.618 1.000 1
## survival_non_censored[74] 74.163 1.000 1
## survival_non_censored[75] 64.873 1.000 1
## survival_non_censored[76] 77.106 1.000 1
## survival_non_censored[77] 56.170 1.000 1
## survival_non_censored[78] 66.585 1.000 1
## survival_non_censored[79] 63.239 1.000 1
## survival_non_censored[80] 70.289 1.000 1
## survival_non_censored[81] 49.834 1.000 1
## survival_non_censored[82] 57.120 1.000 1
## survival_non_censored[83] 71.266 1.000 1
## survival_non_censored[84] 67.559 1.000 1
## survival_non_censored[85] 58.623 1.000 1
## survival_non_censored[86] 60.192 1.000 1
## survival_non_censored[87] 72.474 1.000 1
## survival_non_censored[88] 67.097 1.000 1
## survival_non_censored[89] 57.744 1.000 1
## survival_non_censored[90] 60.504 1.000 1
## survival_non_censored[91] 90.193 1.000 1
## survival_non_censored[92] 63.195 1.000 1
## survival_non_censored[93] 88.632 1.000 1
## survival_non_censored[94] 78.450 1.000 1
## survival_non_censored[95] 46.483 1.000 1
## survival_non_censored[96] 64.405 1.000 1
## survival_non_censored[97] 81.227 1.000 1
## survival_non_censored[98] 93.451 1.000 1
## survival_non_censored[99] 93.328 1.000 1
## survival_non_censored[100] 94.467 1.000 1
## survival_non_censored[101] 69.673 1.000 1
## survival_non_censored[102] 64.277 1.000 1
## survival_non_censored[103] 46.782 1.000 1
## survival_non_censored[104] 35.877 1.000 1
## survival_non_censored[105] 36.235 1.000 1
## survival_non_censored[106] 28.476 1.000 1
## survival_non_censored[107] 39.751 1.000 1
## survival_non_censored[108] 34.856 1.000 1
## survival_non_censored[109] 29.478 1.000 1
## survival_non_censored[110] 60.186 1.000 1
## survival_non_censored[111] 35.337 1.000 1
## survival_non_censored[112] 56.956 1.000 1
## survival_non_censored[113] 38.081 1.000 1
## survival_non_censored[114] 21.769 1.000 1
## survival_non_censored[115] 33.243 1.000 1
## survival_non_censored[116] 22.897 1.000 1
## survival_non_censored[117] 68.624 1.000 1
## survival_non_censored[118] 56.816 1.000 1
## survival_non_censored[119] 87.203 1.000 1
## survival_non_censored[120] 59.800 1.000 1
## survival_non_censored[121] 67.444 1.000 1
## survival_non_censored[122] 30.100 1.000 1
## survival_non_censored[123] 18.710 1.000 1
## survival_non_censored[124] 37.251 1.000 1
## survival_non_censored[125] 61.689 1.000 1
## survival_non_censored[126] 82.998 1.000 1
## survival_non_censored[127] 66.998 1.000 1
## survival_non_censored[128] 44.999 1.000 1
## survival_non_censored[129] 27.907 1.000 1
## survival_non_censored[130] 64.999 1.000 1
## survival_non_censored[131] 23.707 1.000 1
## survival_non_censored[132] 49.002 1.000 1
## survival_non_censored[133] 45.996 1.000 1
## survival_non_censored[134] 9.259 1.000 1
## survival_non_censored[135] 21.506 1.000 1
## survival_non_censored[136] 17.248 1.000 1
## survival_non_censored[137] 18.242 1.000 1
## survival_non_censored[138] 57.843 1.000 1
## survival_non_censored[139] 40.241 1.000 1
## survival_non_censored[140] 51.652 1.000 1
## survival_non_censored[141] 77.352 1.000 1
## survival_non_censored[142] 56.331 1.000 1
## survival_non_censored[143] 27.302 1.000 1
## survival_non_censored[144] 51.526 1.000 1
## survival_non_censored[145] 58.300 1.000 1
## survival_non_censored[146] 63.493 1.000 1
## survival_non_censored[147] 62.667 1.000 1
## survival_non_censored[148] 75.639 1.000 1
## survival_non_censored[149] 45.818 1.000 1
## survival_non_censored[150] 84.203 1.000 1
## survival_non_censored[151] 72.903 1.000 1
## survival_non_censored[152] 21.528 1.000 1
## survival_non_censored[153] 47.220 1.000 1
## survival_non_censored[154] 78.209 1.000 1
## survival_non_censored[155] 34.237 1.000 1
## survival_non_censored[156] 35.318 1.000 1
## survival_non_censored[157] 30.300 1.000 1
## survival_non_censored[158] 21.380 1.000 1
## survival_non_censored[159] 73.248 1.000 1
## survival_non_censored[160] 21.363 1.000 1
## survival_non_censored[161] 69.216 1.000 1
## survival_non_censored[162] 45.936 1.000 1
## survival_non_censored[163] 35.526 1.000 1
## survival_non_censored[164] 59.734 1.000 1
## survival_non_censored[165] 47.316 1.000 1
## survival_non_censored[166] 87.693 1.000 1
## deviance 1433.614 1.010 250
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 11.1 and DIC = 1434.4
## DIC is an estimate of expected predictive error (lower deviance is better).
library(coda)
parameters = c("beta_0", "beta_1", "beta_2", "beta_3", "beta_4", "beta_5",
"beta_6", "beta_7", "beta_8", "beta_9", "r")
res = data.frame(Bayesian_Popes_second_approach$BUGSoutput$sims.matrix)[,parameters]
colnames(res) <- gsub("^.*\\.","", colnames(res) )
# lapply(seq(1,length(res)), function(i) {
# traceplot(as.mcmc(res[,i]), smooth = FALSE, type = "l",
# xlab = "Iterations", ylab = colnames(res)[i],
# main = paste("Traceplot of ", colnames(res)[i]))
# })
tps <- function(var){
ggplot(res, aes_(y=as.name(var), x=seq(1,nrow(res)))) +
geom_line() +
labs(title=paste("Traceplot of ", as.name(var)),
x ="Iterations", y = as.name(var))
}
lapply(names(res), tps)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
##
## [[10]]
##
## [[11]]
library(ggplot2)
res_lag1 = lapply(seq(1:length(res)), function(i) {
lres = lag(res[,i],1)
plot(y=res[,i], x= lres,
xlab = paste0(colnames(res)[i], "-1"),
ylab = paste0(colnames(res)[i]),
main = paste("Lag-1 Scatter Plot of", colnames(res)[i]))
})
lapply(seq(1,length(res)), function(i) {
acf(res[,i], xlab = "Lag", ylab = "ACF",
main = paste("ACF Plot of ", colnames(res)[i]))
})
## [[1]]
##
## Autocorrelations of series 'res[, i]', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.017 0.037 -0.022 0.008 -0.004 0.022 0.022 -0.018 0.001 0.039
## 11 12 13 14 15 16 17 18 19 20 21
## -0.003 -0.004 0.007 -0.018 -0.012 0.003 0.004 -0.001 0.003 -0.027 -0.009
## 22 23 24 25 26 27 28 29 30 31 32
## 0.016 -0.020 -0.004 -0.012 -0.020 0.028 0.010 -0.013 -0.016 0.008 0.010
## 33 34
## 0.000 0.021
##
## [[2]]
##
## Autocorrelations of series 'res[, i]', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.001 0.039 -0.027 0.002 0.016 -0.002 0.006 0.029 -0.014 -0.002
## 11 12 13 14 15 16 17 18 19 20 21
## -0.006 0.020 0.015 -0.005 0.011 -0.021 0.003 -0.016 0.024 -0.003 -0.024
## 22 23 24 25 26 27 28 29 30 31 32
## -0.015 0.005 0.051 -0.022 -0.007 0.027 0.014 0.007 -0.039 -0.019 -0.034
## 33 34
## -0.016 0.006
##
## [[3]]
##
## Autocorrelations of series 'res[, i]', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 -0.005 0.008 0.020 0.003 0.017 0.039 0.000 -0.028 0.011 0.008
## 11 12 13 14 15 16 17 18 19 20 21
## 0.008 0.002 -0.019 -0.002 0.005 0.016 0.016 -0.019 0.021 -0.026 0.008
## 22 23 24 25 26 27 28 29 30 31 32
## -0.001 0.004 0.037 0.007 0.009 0.015 0.004 0.014 0.001 0.015 0.006
## 33 34
## 0.020 0.002
##
## [[4]]
##
## Autocorrelations of series 'res[, i]', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 -0.036 0.029 -0.008 0.016 0.035 0.002 0.009 0.027 -0.019 0.028
## 11 12 13 14 15 16 17 18 19 20 21
## 0.007 0.005 -0.005 0.005 0.002 0.003 -0.023 0.010 -0.004 0.020 0.016
## 22 23 24 25 26 27 28 29 30 31 32
## 0.015 -0.015 0.009 -0.003 -0.034 0.000 0.021 -0.028 -0.025 -0.023 0.017
## 33 34
## 0.007 0.008
##
## [[5]]
##
## Autocorrelations of series 'res[, i]', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.024 -0.001 0.010 0.005 -0.030 0.005 0.016 0.007 -0.010 0.003
## 11 12 13 14 15 16 17 18 19 20 21
## -0.006 -0.032 -0.009 0.009 -0.008 -0.020 -0.016 -0.027 -0.017 -0.006 -0.004
## 22 23 24 25 26 27 28 29 30 31 32
## -0.011 0.004 0.039 0.031 -0.022 -0.018 -0.002 0.011 0.002 0.005 -0.004
## 33 34
## -0.029 0.026
##
## [[6]]
##
## Autocorrelations of series 'res[, i]', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.002 -0.020 -0.026 0.005 -0.005 0.009 0.014 0.016 0.007 0.016
## 11 12 13 14 15 16 17 18 19 20 21
## -0.017 -0.028 0.016 0.002 0.001 -0.030 -0.002 -0.019 0.000 -0.006 0.012
## 22 23 24 25 26 27 28 29 30 31 32
## 0.009 -0.038 0.004 0.008 -0.019 0.007 -0.035 0.024 -0.012 -0.010 -0.002
## 33 34
## 0.025 0.021
##
## [[7]]
##
## Autocorrelations of series 'res[, i]', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 -0.026 0.012 0.008 0.023 -0.010 0.014 -0.009 -0.030 -0.011 0.007
## 11 12 13 14 15 16 17 18 19 20 21
## 0.015 0.001 -0.021 0.028 0.007 -0.010 0.013 -0.013 0.034 -0.034 0.011
## 22 23 24 25 26 27 28 29 30 31 32
## 0.012 0.016 0.026 0.010 -0.010 0.000 0.031 0.024 -0.004 0.015 0.002
## 33 34
## 0.016 0.038
##
## [[8]]
##
## Autocorrelations of series 'res[, i]', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 -0.015 0.013 -0.003 -0.008 -0.014 -0.037 -0.013 0.005 0.021 0.001
## 11 12 13 14 15 16 17 18 19 20 21
## -0.008 0.043 0.009 0.019 -0.017 0.023 -0.021 -0.002 -0.013 -0.001 -0.035
## 22 23 24 25 26 27 28 29 30 31 32
## -0.027 0.012 0.034 0.020 0.010 -0.006 0.012 0.002 -0.008 0.009 -0.022
## 33 34
## -0.013 0.011
##
## [[9]]
##
## Autocorrelations of series 'res[, i]', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 -0.023 0.039 -0.012 -0.005 -0.022 -0.038 -0.003 0.022 -0.026 -0.015
## 11 12 13 14 15 16 17 18 19 20 21
## 0.001 0.008 0.012 0.015 0.003 0.015 -0.031 -0.003 0.018 0.006 0.021
## 22 23 24 25 26 27 28 29 30 31 32
## 0.000 -0.036 0.018 0.008 -0.025 0.031 0.018 0.006 -0.005 -0.033 0.007
## 33 34
## 0.007 -0.015
##
## [[10]]
##
## Autocorrelations of series 'res[, i]', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.022 0.022 0.006 -0.001 0.019 0.002 0.035 0.016 -0.012 -0.006
## 11 12 13 14 15 16 17 18 19 20 21
## -0.016 0.021 -0.006 -0.016 0.009 -0.003 -0.011 0.012 0.023 0.006 -0.053
## 22 23 24 25 26 27 28 29 30 31 32
## 0.028 -0.002 0.022 -0.029 -0.017 0.011 -0.016 0.012 -0.008 -0.043 -0.023
## 33 34
## -0.026 0.014
##
## [[11]]
##
## Autocorrelations of series 'res[, i]', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10
## 1.000 0.014 0.011 -0.013 0.008 0.015 0.003 0.013 -0.014 0.002 0.027
## 11 12 13 14 15 16 17 18 19 20 21
## 0.007 -0.013 0.014 -0.009 -0.010 0.009 0.012 0.028 0.005 0.009 -0.020
## 22 23 24 25 26 27 28 29 30 31 32
## -0.011 -0.012 -0.019 -0.018 -0.020 0.024 -0.006 -0.003 0.006 0.012 0.003
## 33 34
## -0.016 0.002
people = c("age_Francis_predictive", "age_Obama_predictive", "age_Dalai_predictive",
"age_Naruhito_predictive")
four_pred = data.frame(Bayesian_Popes_second_approach$BUGSoutput$sims.matrix)[,people]
print(paste("P(Dalai Lifespan > Francis Lifespan)=", mean(four_pred$age_Dalai_predictive>four_pred$age_Francis_predictive)))
## [1] "P(Dalai Lifespan > Francis Lifespan)= 0.383666666666667"
print(paste("P(Obama Lifespan > Naruhito Lifespan)=", mean(four_pred$age_Obama_predictive>four_pred$age_Naruhito_predictive)))
## [1] "P(Obama Lifespan > Naruhito Lifespan)= 0.501"
ppd_plot <- function(i) {
name = gsub(".*[_]([^.]+)[_].*", "\\1", colnames(four_pred)[i])
minval = round(min(four_pred[,i]), 2)
medval = round(median(four_pred[,i]), 2)
val95 = round(quantile(four_pred[,i], 0.95),2)
maxval = round(max(four_pred[,i]), 2)
ggplot(four_pred, aes(x=four_pred[,i])) + geom_histogram(aes(y=..density..),
binwidth = 0.8,
boundary = minval) +
geom_vline(xintercept = minval, colour = "black")+
geom_vline(xintercept = medval, colour = "black") +
geom_vline(xintercept = val95, colour = "black")+
geom_hline(yintercept = 0, colour = "black") +
xlab("Lifespan (years)") +
ylab("Predictive Probability density") +
ggtitle(paste("Predictive Posterior Distribution for Lifespan of", name)) +
ylim(c(0, 0.07))+
annotate("text", x = minval+0.15*(medval-minval), y = 0.06, label = paste("Current age", minval, "years", sep = "\n"), size = 3)+
annotate("text", x = (medval+minval)/2, y = 0.05, label = "50%", size = 3)+
annotate("text", x = medval+0.15*(val95-medval), y = 0.06, label = paste("Median", medval,"years", sep = "\n"), size = 3)+
annotate("text", x = (medval+val95)/2, y = 0.05, label = "45%", size = 3)+
annotate("text", x = val95+2.2, y = 0.06, label = paste("95%", "quantile", val95, "years", sep = "\n"), size = 3)+
annotate("text", x = val95+3, y = 0.05, label = "5%", size = 3)
}
lapply(seq(1,4), ppd_plot)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
francis_dalai = data.frame(plifespans = c(four_pred$age_Francis_predictive, four_pred$age_Dalai_predictive), person = c(rep("Francis", nrow(four_pred)), rep("Dalai", nrow(four_pred))))
mins = francis_dalai %>%
group_by(person) %>%
summarise(min_vals = round(min(plifespans), 2))
## `summarise()` ungrouping output (override with `.groups` argument)
meds = francis_dalai %>%
group_by(person) %>%
summarise(med_vals = round(median(plifespans), 2))
## `summarise()` ungrouping output (override with `.groups` argument)
val95 = francis_dalai %>%
group_by(person) %>%
summarise(vals_95 = round(quantile(plifespans, 0.95), 2))
## `summarise()` ungrouping output (override with `.groups` argument)
francis_txt = paste(paste("Current age:", mins$min_vals[mins$person=="Francis"], "yrs"),
paste("Median:", meds$med_vals[meds$person=="Francis"],"yrs"),
paste("95% quantile:", val95$vals_95[val95$person=="Francis"], "yrs"),
sep = "\n")
dalai_txt = paste(paste("Current age:", mins$min_vals[mins$person=="Dalai"], "yrs"),
paste("Median:", meds$med_vals[meds$person=="Dalai"],"yrs"),
paste("95% quantile:", val95$vals_95[val95$person=="Dalai"], "yrs"),
sep = "\n")
plot_txt <- data.frame(
label = c(francis_txt, dalai_txt),
person = c("Francis", "Dalai")
)
ggplot(francis_dalai, aes(x=plifespans)) +
geom_histogram(aes(y=..density..), binwidth = 0.8,boundary = round(min(francis_dalai$plifespans), 2)) +
facet_wrap(~person) +
geom_vline(data=mins, aes(xintercept=min_vals))+
geom_vline(data=meds, aes(xintercept=med_vals))+
geom_vline(data=val95, aes(xintercept=vals_95))+
geom_hline(yintercept = 0, colour = "black") +
xlab("Lifespan (yrs)") + ylab("Predictive Probability density") +
ggtitle(paste("Posterior Predictive Lifespans (Pope Francis and 14th Dalai Lama)")) +
ylim(c(0, 0.05)) +
geom_text(data = plot_txt, mapping = aes(x = 92.8, y = 0.048, label = label), size = 2.8)
obama_naru = data.frame(plifespans = c(four_pred$age_Obama_predictive, four_pred$age_Naruhito_predictive), person = c(rep("Obama", nrow(four_pred)), rep("Naruhito", nrow(four_pred))))
mins = obama_naru %>%
group_by(person) %>%
summarise(min_vals = round(min(plifespans), 2))
## `summarise()` ungrouping output (override with `.groups` argument)
meds = obama_naru %>%
group_by(person) %>%
summarise(med_vals = round(median(plifespans), 2))
## `summarise()` ungrouping output (override with `.groups` argument)
val95 = obama_naru %>%
group_by(person) %>%
summarise(vals_95 = round(quantile(plifespans, 0.95), 2))
## `summarise()` ungrouping output (override with `.groups` argument)
obama_txt = paste(paste("Current age:", mins$min_vals[mins$person=="Obama"], "yrs"),
paste("Median:", meds$med_vals[meds$person=="Obama"],"yrs"),
paste("95% quantile:", val95$vals_95[val95$person=="Obama"], "yrs"),
sep = "\n")
naru_txt = paste(paste("Current age:", mins$min_vals[mins$person=="Naruhito"], "yrs"),
paste("Median:", meds$med_vals[meds$person=="Naruhito"],"yrs"),
paste("95% quantile:", val95$vals_95[val95$person=="Naruhito"], "yrs"),
sep = "\n")
plot_txt <- data.frame(
label = c(obama_txt, naru_txt),
person = c("Obama", "Naruhito")
)
ggplot(obama_naru, aes(x=plifespans)) +
geom_histogram(aes(y=..density..), binwidth = 0.8,boundary = round(min(obama_naru$plifespans), 2)) +
facet_wrap(~person) +
geom_vline(data=mins, aes(xintercept=min_vals))+
geom_vline(data=meds, aes(xintercept=med_vals))+
geom_vline(data=val95, aes(xintercept=vals_95))+
geom_hline(yintercept = 0, colour = "black") +
xlab("Lifespan (yrs)") + ylab("Predictive Probability density") +
ggtitle(paste("Posterior Predictive Lifespans (President Obama and Emperor Naruhito)")) +
ylim(c(0, 0.05)) +
geom_text(data = plot_txt, mapping = aes(x = 80, y = 0.04, label = label), size = 3)
vars = data.frame(Bayesian_Popes_second_approach$BUGSoutput$root.short)
output <- head(data.frame(Bayesian_Popes_second_approach$BUGSoutput$summary), 16)
output <- output %>%
mutate("Variable" = head(vars,16),
"Mean" = output$mean,
"Standard Deviation" = output$sd,
"2.5% Quantile" = output$X2.5.,
"Median" = output$X50.,
"97.5% Quantile" = output$X97.5.
) %>%
select("Variable", "Mean", "Standard Deviation", "2.5% Quantile", "Median", "97.5% Quantile")
coef <- output %>%
slice(5:16)
kable(coef)
| Variable | Mean | Standard Deviation | 2.5% Quantile | Median | 97.5% Quantile |
|---|---|---|---|---|---|
| beta_0 | -22.7763077 | 1.5962011 | -25.9407203 | -22.7453356 | -19.5873392 |
| beta_1 | -0.0016375 | 0.0007812 | -0.0031670 | -0.0016413 | -0.0001422 |
| beta_2 | 0.6008876 | 0.4402710 | -0.2822430 | 0.6097154 | 1.4157974 |
| beta_3 | 1.3236079 | 0.2677447 | 0.7903557 | 1.3239796 | 1.8480258 |
| beta_4 | 1.6942849 | 0.3354591 | 1.0192361 | 1.7074271 | 2.3307992 |
| beta_5 | 0.7160662 | 0.2440568 | 0.2412665 | 0.7124606 | 1.2012485 |
| beta_6 | -0.0020120 | 0.0024374 | -0.0067747 | -0.0020312 | 0.0026784 |
| beta_7 | 0.0006767 | 0.0013981 | -0.0020156 | 0.0006623 | 0.0034650 |
| beta_8 | 0.0068935 | 0.0017716 | 0.0036422 | 0.0068926 | 0.0103842 |
| beta_9 | 0.0002039 | 0.0012350 | -0.0021660 | 0.0001918 | 0.0025679 |
| deviance | 1423.3093231 | 4.7248506 | 1416.0482075 | 1422.6535327 | 1433.6141160 |
| r | 4.2402583 | 0.2647113 | 3.7428743 | 4.2410108 | 4.7623073 |
beta_0_est <- coef[1,2]
beta_1_est <- coef[2,2]
beta_2_est <- coef[3,2]
beta_3_est <- coef[4,2]
beta_4_est <- coef[5,2]
beta_5_est <- coef[6,2]
beta_6_est <- coef[7,2]
beta_7_est <- coef[8,2]
beta_8_est <- coef[9,2]
beta_9_est <- coef[10,2]
no_living <- leaders_data %>%
filter(Censored == 0)
estimates <- no_living %>%
mutate(Birth.Year.Cent = Birth.Year - mean(leaders_nopred$Birth.Year)) %>%
mutate(ChinaEmp =case_when(Type == "ChinaEmp" ~ 1, TRUE ~ 0),
DalaiLama =case_when(Type == "DalaiLama" ~ 1, TRUE ~ 0),
JapanEmp =case_when(Type == "JapanEmp" ~ 1, TRUE ~ 0),
Pope =case_when(Type == "Pope" ~ 1, TRUE ~ 0),
UsPres =case_when(Type == "UsPres" ~ 1, TRUE ~ 0)) %>%
mutate(BYUsPres = Birth.Year.Cent*UsPres) %>%
mutate(BYChinaEmp = Birth.Year.Cent*ChinaEmp) %>%
mutate(BYDalaiLama = Birth.Year.Cent*DalaiLama) %>%
mutate(BYJapanEmp = Birth.Year.Cent*JapanEmp) %>%
mutate(intercept = 1) %>%
mutate(predicted = exp(beta_0_est + beta_1_est*(Birth.Year.Cent) + beta_2_est*(UsPres) + beta_3_est*(ChinaEmp) + beta_4_est*(DalaiLama) + beta_5_est*(JapanEmp) + beta_6_est*(UsPres)*(Birth.Year.Cent) + beta_7_est*(ChinaEmp)*(Birth.Year.Cent) + beta_8_est*(DalaiLama)*(Birth.Year.Cent) + beta_9_est*(JapanEmp)*(Birth.Year.Cent)))
# check using matrices -- Ts in the 200s
data.mat = as.matrix(estimates[,c("intercept", "Birth.Year.Cent", "UsPres", "ChinaEmp", "DalaiLama", "JapanEmp","BYUsPres", "BYChinaEmp", "BYDalaiLama", "BYJapanEmp")], nrow = nrow(estimates), byrow = TRUE)
r = Bayesian_Popes_second_approach$BUGSoutput$summary[,1][16]
betas = Bayesian_Popes_second_approach$BUGSoutput$summary[,1][5:14]
beta.mat = as.matrix(betas, nrow = p, byrow = TRUE)
logmus = data.mat %*% beta.mat
min.age = min(leaders_nopred$Age.Event)
install.packages("TruncatedDistributions", repos="http://R-Forge.R-project.org")
## Installing package into '/home/rstudio-user/R/x86_64-pc-linux-gnu-library/4.0'
## (as 'lib' is unspecified)
library(TruncatedDistributions)
Ts = rtweibull(166, r, (1/exp(logmus))^(1/r), min.age, 120)
# dweibull takes a diff parameterization of Weibull
# please check the parameterization see ?dweibull
non_censored_prediction <- filter(leaders_data, Censored == 0)
K <- 10
ysims <- matrix(nrow = nrow(res), ncol = K)
for(k in 1:K){
for(i in 1:nrow(non_censored_prediction)){
w <- sample(nrow(res), 1)
samp <- res[w,]
val <- non_censored_prediction[i, ]
beta_temp <- samp$beta_0 + samp$beta_1*val$Birth.Year + samp$beta_2*as.integer(val$TypeUsPres) + samp$beta_3*as.integer(val$TypeChinaEmp) + samp$beta_4*as.integer(val$TypeDalaiLama) + samp$beta_5*as.integer(val$TypeJapanEmp) + samp$beta_6*val$Birth.Year*as.integer(val$TypeUsPres) + samp$beta_7*val$Birth.Year*as.integer(val$TypeChinaEmp) + samp$beta_8*val$Birth.Year*as.integer(val$TypeDalaiLama) + samp$beta_9*val$Birth.Year*as.integer(val$TypeJapanEmp)
mu <- exp(beta_temp)
t <- dweibull(samp$r, mu)
ysims[i, k] <- t
}
}
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
## Warning in dweibull(samp$r, mu): NaNs produced
r <- sample(nrow(non_censored_prediction), 1)
hist(ysims[r,])
abline(v=non_censored_prediction$Age.Event[r], col = "red")
print(mean(ysims<0))
## [1] NA
print(mean(non_censored_prediction$Age.Event <0))
## [1] 0
#lg <- lexis_grid(1960, 2000, 30, 80, delta=5)
#leaders_data$Initial.Age <- year(leaders_data$Event.Date) - leaders_data$Birth.Year
#subset <- filter(leaders_data, Censored == 0)
#subset <- subset[0:5,]
#ggplot() + geom_polygon(data = subset, mapping = aes(x = Birth.Year, y = Age.Event, group = Name)) + geom_point(data = subset, aes(x = Birth.Year, y = Birth.Year + Age.Event)) + geom_polygon(data = subset, mapping = aes(x = Birth.Year, y = ))
subset <- filter(leaders_data, Censored == 0)
ggplot(data = subset, aes(x = Birth.Year, Age.Event)) + geom_point(color="pink") + ggtitle("Distribution of Leaders' Ages") + labs(x = "Birth Year", y = "Age")